
#################################################################################
# Fig. 5. Effect of voting threshold on winning coalition size by membership size
#################################################################################

# Programme:  cbld-analysis02-fig05.r
# Project:		Council coalition building
# Author:		  Frank Haege, Department of Politics and Administration, University of Limerick
# Contact:		frank.haege@ul.ie

# Description
#############
# This R-script takes the simulated consensus rates data exported from Stata
# and fits a non-linear quantile regression for each membership size
# The results are imported into Stata for further analysis through the same
# do-file ('cbld-analysis02-fig05.do')


# Remove all objects
rm(list = ls(all = TRUE))

# Load libraries
library(foreign)
library(quantreg)

# Set working directory
setwd("C:\\")

# Load stata dataset
exp01 = read.dta("Experiments\\Experiment01\\cbld-analysis02-fig05.dta", convert.dates = TRUE)

# Estimate non-linear quantile regression for each membership size
for (i in c(6, 9, 12, 15, 18, 21, 24, 27, 30)) { 
  # Generate indicator variable for relevant sub-sample
  selms = exp01$nostates == i
  # Extract the relevant subsample and name the new dataframe according to the loop counter
  dta = assign(paste('ms', i, sep=''), exp01[selms, ])
  # Attach the current dataframe
  attach(dta)
  # Fit the simple product function through a non-linear quantile regression
  fit1lq = nlrq(normconsensus ~ 0 + I(1-(1-normvthresh)^exponent), tau = 0.25, start = list(exponent = 0.25), data = dta)
  fit1m = nlrq(normconsensus ~ 0 + I(1-(1-normvthresh)^exponent), tau = 0.5, start = list(exponent = 0.5), data = dta)
  fit1uq = nlrq(normconsensus ~ 0 + I(1-(1-normvthresh)^exponent), tau = 0.75, start = list(exponent = 0.75), data = dta)
  # Extract the coefficient estimate and save it as a new variable in the dataframe
  dta$cfit1lq = coef(fit1lq)
  dta$cfit1m = coef(fit1m)
  dta$cfit1uq = coef(fit1uq)    
  # Extract the residuals and save them in a new variable in the dataframe
  dta$rfit1lq = residuals(fit1lq)
  dta$rfit1m = residuals(fit1m)
  dta$rfit1uq = residuals(fit1uq)
  # Extract the predicted values and save them in a new variable in the dataframe
  dta$pfit1lq = predict(fit1lq)
  dta$pfit1m = predict(fit1m)
  dta$pfit1uq = predict(fit1uq)
  # Fit the product function with intercept
  fit2lq = nlrq(normconsensus ~ 0 + intercept + I((1-intercept)*(1-(1-normvthresh)^exponent)), tau = 0.25, start = list(intercept = 0.3, exponent = 10), data = dta)
  fit2m = nlrq(normconsensus ~ 0 + intercept + I((1-intercept)*(1-(1-normvthresh)^exponent)), tau = 0.5, start = list(intercept = 0.3, exponent = 10), data = dta)
  fit2uq = nlrq(normconsensus ~ 0 + intercept + I((1-intercept)*(1-(1-normvthresh)^exponent)), tau = 0.75, start = list(intercept = 0.3, exponent = 10), data = dta)
  # Extract the intercept estimate and save it as a new variable in the dataframe
  dta$ifit2lq = coef(fit2lq)[1]
  dta$ifit2m = coef(fit2m)[1]
  dta$ifit2uq = coef(fit2uq)[1]
  # Extract the coefficient estimate and save it as a new variable in the dataframe
  dta$cfit2lq = coef(fit2lq)[2]
  dta$cfit2m = coef(fit2m)[2]
  dta$cfit2uq = coef(fit2uq)[2]  
  # Extract the residuals and save them in a new variable in the dataframe
  dta$rfit2lq = residuals(fit2lq)
  dta$rfit2m = residuals(fit2m)
  dta$rfit2uq = residuals(fit2uq)
  # Extract the predicted values and save them in a new variable in the dataframe
  dta$pfit2lq = predict(fit2lq)
  dta$pfit2m = predict(fit2m)
  dta$pfit2uq = predict(fit2uq)
  # Update the original dataframe with the new variables
  assign(paste('ms', i, sep=''), dta)
  # Detach the current dataframe
  detach(dta)
}

# Stack datasets with predicted values for voting threshold effect 
ms = rbind(ms6, ms9, ms12, ms15, ms18, ms21, ms24, ms27, ms30)
rm(dta)

# Estimate linear quantile regression for each voting threshold
for (i in c(51, 55, 60, 65, 70, 75, 80, 85, 90)) { 
  # Generate indicator variable for relevant sub-sample
  selvt = ms$vthresh == i
  # Extract the relevant subsample and name the new dataframe according to the loop counter
  dta = assign(paste('vt', i, sep=''), ms[selvt, ])
  # Attach the current dataframe
  attach(dta)
  # Fit the linear quantile regression
  fit3lq = rq(consensus ~ nostates, tau = 0.25, data = dta)
  fit3m = rq(consensus ~ nostates, tau = 0.5, data = dta)
  fit3uq = rq(consensus ~ nostates, tau = 0.75, data = dta)
  # Extract the intercept estimate and save it as a new variable in the dataframe
  dta$ifit3lq = coef(fit3lq)[1]
  dta$ifit3m = coef(fit3m)[1]
  dta$ifit3uq = coef(fit3uq)[1]  
  # Extract the coefficient estimate and save it as a new variable in the dataframe
  dta$cfit3lq = coef(fit3lq)[2]
  dta$cfit3m = coef(fit3m)[2]
  dta$cfit3uq = coef(fit3uq)[2]    
  # Extract the residuals and save them in a new variable in the dataframe
  dta$rfit3lq = residuals(fit3lq)
  dta$rfit3m = residuals(fit3m)
  dta$rfit3uq = residuals(fit3uq)
  # Extract the predicted values and save them in a new variable in the dataframe
  dta$pfit3lq = predict(fit3lq)
  dta$pfit3m = predict(fit3m)
  dta$pfit3uq = predict(fit3uq)
  # Update the original dataframe with the new variables
  assign(paste('vt', i, sep=''), dta)
  # Detach the current dataframe
  detach(dta)
}

# Stack datasets with predicted values for voting threshold effect 
results = rbind(vt51, vt55, vt60, vt65, vt70, vt75, vt80, vt85, vt90)

# Save estimation results in stata format
write.dta(results, "Experiments\\Experiment01\\cbld-analysis02-fig05-r.dta") 